Temnos2021Analysis
Phoebe Koenig
2023-03-09
Loading
Packages
library(lme4)## Loading required package: Matrixlibrary(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmListlibrary(car)## Loading required package: carDatalibrary(ggplot2)
library(plyr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following object is masked from 'package:car':
##
## recode## The following object is masked from 'package:nlme':
##
## collapse## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBalibrary(performance)
library(patchwork)##
## Attaching package: 'patchwork'## The following object is masked from 'package:MASS':
##
## arealibrary(emmeans)
library(MuMIn)
library(rsq)Data
Import dataset and remove ants without measured venom
datastore <- read.csv("~/Temnos2021/Koenig_TemnosVenom_2021.csv")
data<-na.omit(datastore)Variables
List:
Dataset has dissected workers and queens. It also has some pupae. Here are the variables in the dataset:
Date.Collected.D.M : (DD-Month-YYYY) Date on which colony was collected. Year should always be 2021.
Date.Dissected.D.M : (DD-Month-YYYY) Date on which ant was dissected to measure venom sac. Year should be 2021 or 2022.
Date.Collected.Year : (YYYY) Year in which colony was collected. Should be 2021 for every colony.
Colony : The colony number identifier. Technically, this is a nest, not a colony.
Adult.Count : The number of adults in the nest on the day it was censused
Queen.Count : The number of dealate queens (wingless) in the nest on the day it was censused
Larvae.Count : The number of larvae in the nest on the day it was censused.
Pupae.Count : The number of pupae in the nest on the day it was censused Males.Count : The number of males in the nest on the day it was censused
Winged.Queen.Count : The number of alate queens (winged) in the nest on the day it was censused
Egg.Count : The number of eggs in the nest on the day it was censused. This number could be an underestimate as eggs are small and hard to count
Worker.Queen : Whether the dissected ant was a worker or a queen
Behavior.Code : If Worker.Queen==“Worker”, this field should be either Nurse or Forager. If Worker.Queen==“Queen”, this field should be either Mated or Unmated. Mated means the queen was counted in the Queen.Count and was dealate. Unmated means the queen was either counted in Winged.Queen.Count or hatched while the colony was in the lab, and was alate. For both workers and queens, we also have the option pupa, which is where I dissected a pupa. There are not very many of them and I exclude these data for analyses.
Callow : (Binary) was the individual dissected callow? 1 if yes, 0 if not
Webers.Length.mm : Weber’s length in mm
Venom.Sac.Length.mm : Length of venom sac in mm
Venom.Sac.Width.mm : Width of venom sac in mm
Ant : Ant number (identifier)
Sac.Collected : (binary) whether or not the sac was successfully harvested to store in freezer for chemical analysis? 1 if yes, 0 if no
Add and mutate:
In this chunk, I import the summary datasheet (nest level data) and add a few new variables to the dissection dataset that can be calculated from the initial variables. Almost all info included in the summary dataset is included in the original dataset, but I import the summary dataset in order to scale the nest census variables, because I want them scaled by the average colony size. I also make sure the date is in the correct format. The most important variable created here is Venom.Volume, which is calculated with the formula Venom.Volume \(=\frac{\pi}{6} \times(L \times W^2)\) for each venom sac.
I use Date to create a new variable, Collection.Day where day 1 is the first day I collected ants in the field.
I subset the dataset to create separate datasets of workers and queens.
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, uniondata$Date.Collected.D.M <- as.Date((data$Date.Collected.D.M), format = "%d-%b-%y")
data$Date.Dissected.D.M <- as.Date((data$Date.Dissected.D.M), format = "%d-%b-%y")
data$Callow<-as.logical(data$Callow)
summarydata<- read.csv("~/Temnos2021/Koenig_Temnos_2021_nest_census_data.csv")
summarydata$Date.Collected <- as.Date(as.character(summarydata$Date.Collected), format = "%d- %b")
year(summarydata$Date.Collected) <- 2021
#R auto assigned current year, but these nests were all collected in 2021, so I reassign year
#In R, first day is automatically Jan 1, 1970. therefore we change to day of collection to make results more interpretable. To do this I first convert Collection date to a numeric, where day 1 is the first day I collected ant nests (since first day in R is 1970 I convert to a numeric and then subtract every day before that, which is 18750 days)
data$Collection.Day<-as.numeric(data$Date.Collected.D.M)-18750
#Scale variables to make results more interpretable. Since there are multiple rows of data in our main dataset for the same nest, we don't want nests that had more workers to contribute more to the mean for adult count or larva count. We want only one contribution to the mean per nest. To do this I pulled in the summarydata dataset which is all colonies collected over the season (some of which never led to successful dissections of workers). I then subset the dataset to make sure I'm only calculating means from nests that did see successful dissections
# subset row by colony numbers that are included in data df
summarydatasubset<-subset(summarydata, summarydata$Colony %in% data$Colony)
#calculate means and SDs for Adult, pupa, and larva count. Define scaled variables in data df
meanadults<-mean(summarydatasubset$Adult.Count)
SDadults<-sd(summarydatasubset$Adult.Count)
data$Adult.Count.Scaled<-(data$Adult.Count-meanadults)/SDadults
meanpupae<-mean(summarydatasubset$Pupae.Count)
SDpupae<-sd(summarydatasubset$Pupae.Count)
data$Pupae.Count.Scaled<-(data$Pupae.Count-meanpupae)/SDpupae
meanlarvae<-mean(summarydatasubset$Larvae.Count)
SDlarvae<-sd(summarydatasubset$Larvae.Count)
data$Larvae.Count.Scaled<-(data$Larvae.Count-meanlarvae)/SDlarvae
data$Location.Collected<-as.factor(data$Location.Collected)
#Create a variable for venom volume, which is an ellipsoid using venom sac length and width. Since both length and width are measured in mm, this volume is in mm^3, which is 1:1 with microliters. For the purposes of graphing and modeling this problem, we will multiply by 1000 to convert to nanoliters
data$Venom.Volume<-(pi/6)*data$Venom.Sac.Length.mm*data$Venom.Sac.Width.mm^2*1000
#Standardize venom by body size
data$standardizedvenom<-((data$Venom.Volume)/(data$Webers.Length.mm))
#Create a dataset for just workers. At this point these datasets still contain pupae, which must be removed for analyses
workerdata<-subset(data,Worker.Queen=="Worker")
queendata<-subset(data,Worker.Queen=="Queen")Dataset size
nrow(subset(workerdata,Behavior.Code=="Pupa"))## [1] 15#Dissected 15 worker pupae
nrow(subset(queendata,Behavior.Code=="Pupa"))## [1] 9#Dissected 9 queen pupae
length(unique(workerdata[["Colony"]]))## [1] 155length(unique(queendata[["Colony"]]))## [1] 104#Remove pupae from datasets!!!
workerdata<-subset(workerdata,Behavior.Code!="Pupa")
queendata<-subset(queendata,Behavior.Code!="Pupa")
nrow(workerdata)## [1] 2055#Successfully dissected 2054 adult worker
nrow(subset(workerdata,Behavior.Code=="Nurse"))## [1] 1394length(unique(subset(workerdata,Behavior.Code=="Nurse")$Colony))## [1] 153nrow(subset(workerdata,Behavior.Code=="Forager"))## [1] 661length(unique(subset(workerdata,Behavior.Code=="Forager")$Colony))## [1] 142#1394 of these were nurses from 153 colonies, 661 of them foragers from 142 colonies
nrow(queendata)## [1] 199#Successfully dissected 199 queens
nrow(subset(queendata,Behavior.Code=="Mated"))## [1] 155nrow(subset(queendata,Behavior.Code=="Unmated"))## [1] 44#155 of them were dealate queens, 44 were alate queens. They are coded in dataset as mated or unmated but I did not dissect spermathecas, just looked at whether or not they had wings
nrow(subset(workerdata,Callow==T))## [1] 16#16
length(unique(subset(workerdata,Callow==T)$Colony))## [1] 8# from 8 colonies
nrow(subset(workerdata,Callow==F))## [1] 2039length(unique(subset(workerdata,Callow==F)$Colony))## [1] 155#2039 from 155 colonies
mean(workerdata$Venom.Volume)## [1] 1.182956sd(workerdata$Venom.Volume)## [1] 0.5586431#What is the average number of nests per collection day?
uniquedays<-unique(workerdata$Collection.Day)
length(uniquedays)## [1] 40uniquenests<-unique(workerdata$Colony)
length(uniquenests)## [1] 155length(uniquenests)/length(uniquedays)## [1] 3.875#3.9 nests on average per day, 40 collection days
queened<-subset(workerdata,Queen.Count>0)
notqueened<-subset(workerdata,Queen.Count==0)
length(unique(queened$Colony))## [1] 112length(unique(notqueened$Colony))## [1] 43Visualizing worker data
Behavior Code
At this point, we need to remove pupae from dataset.
workerdata<-subset(workerdata,Behavior.Code!="Pupa")
queendata<-subset(queendata,Behavior.Code!="Pupa")Do nurses and foragers have the same amount of venom? Do workers in queenright colonies have the same amount of venom as workers in queenless colonies? Do callow workers have the same amount of venom as melanized workers? Graph it.
venombybehavior<-ggplot(workerdata,aes(Behavior.Code,Venom.Volume,color=Behavior.Code))+
geom_violin()+
geom_boxplot(width=0.1)+
xlab("Behavioral Group")+ylab("Venom Volume (nL)")+theme_bw(base_size = 22)+
theme(legend.position="none")+
scale_colour_discrete("")
venombyqueenright<-ggplot(workerdata,aes(as.factor(Queen.Count>0),Venom.Volume, color=as.factor(Queen.Count>0)))+
geom_violin()+
geom_boxplot(width=0.1)+
xlab("Number of Dealate Queens")+ylab("Venom Volume (nL)")+theme_bw(base_size = 22)+
theme(legend.position="none")+
scale_x_discrete(labels = c('0','At least 1'))
callowplot = ggplot(data=workerdata,aes(x=Callow, y=Venom.Volume, color=Callow)) + geom_violin() +geom_boxplot(width=0.1)+
xlab("Callow")+ylab("Venom Volume (nL)")+theme_bw(base_size = 22)+theme(legend.position = "none")+ scale_x_discrete(labels=c("No", "Yes"))
remove_y <- theme(
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank()
)
p=list(venombybehavior,callowplot+remove_y,venombyqueenright+remove_y)
wrap_plots(p, nrow = 1) 
ggsave('workers2021.png',plot=wrap_plots(p, nrow = 1),dpi=900,width =39, height = 15, units = "cm",device='png')Worker venom volume represented graphically
#now the initial model, which looks for an interaction effect between behavior and pupae count, and between behavior and adult count.
#exclude dissected pupae from dataset (we already did above, but no problem to have code 2x in case do not run above chunk)
workerdata<-subset(workerdata, Behavior.Code!="Pupa")
hist(workerdata$Venom.Volume)
#A bit right skewed
callowplot = ggplot(data=workerdata,aes(x=Callow, y=standardizedvenom, color=Callow)) + geom_violin() +geom_boxplot(width=0.1)+
xlab("Callow")+ylab("Standardized Venom Volume")+theme_bw(base_size = 18)+theme(legend.position = "none")+ scale_x_discrete(labels=c("No", "Yes"))
callowplot
venombysize = ggplot(data=workerdata,aes(x=Webers.Length.mm, y=Venom.Volume)) + geom_point() +
xlab("Weber's Length (mm)")+ylab("Venom Volume")+theme_bw(base_size = 22)+geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE)
venombysize## `geom_smooth()` using formula = 'y ~ x'
venombydate = ggplot(data=workerdata, aes(Date.Collected.D.M,Venom.Volume)) +
geom_point() +
scale_x_date(date_breaks = "months" , date_labels = "%b")+
geom_smooth(method=lm, se=TRUE ) +
xlab("Date of Nest Collection ")+
ylab("Venom Volume (nL)")+
theme_bw(base_size = 22)
venombydate## `geom_smooth()` using formula = 'y ~ x'
ggsave('venombydate2021revised.png',plot=venombydate,dpi=900,width =25, height = 15, units = "cm",device='png')## `geom_smooth()` using formula = 'y ~ x'venombypupae = ggplot(data=workerdata, aes(Pupae.Count,standardizedvenom)) +
geom_point() +
geom_smooth(method=lm, color="red", fill="gray", se=TRUE) +
xlab("Number of Pupae in Nest")+
ylab("Standardized Venom Volume")+
theme_bw(base_size = 22)
venombypupae## `geom_smooth()` using formula = 'y ~ x'
venombyadultinteraction=ggplot(data=workerdata, aes(Adult.Count,Venom.Volume, color=Behavior.Code)) +
geom_point() +
geom_smooth(method=lm, fill="gray", se=TRUE) +
xlab("Number of Workers in Nest")+
ylab("Venom Volume (nL)")+
theme_bw(base_size = 22)+
scale_colour_discrete("Behavior")
venombyadultinteraction## `geom_smooth()` using formula = 'y ~ x'
ggsave('adultinteraction2021revised.png',plot=venombyadultinteraction,dpi=900,width =30, height = 15, units = "cm",device='png')## `geom_smooth()` using formula = 'y ~ x'#Define ratio of pupae to adults
workerdata$Pupae.Adult.Ratio=workerdata$Pupae.Count/workerdata$Adult.Count
venombyratio=ggplot(data=workerdata, aes(Pupae.Adult.Ratio,Venom.Volume, color=Behavior.Code)) +
geom_point() +
geom_smooth(method=lm, fill="gray", se=TRUE) +
xlab("Ratio of Pupae to Adults")+
ylab("Venom Volume (nL)")+
theme_bw(base_size = 22)+
scale_colour_discrete("Behavior")
venombyratio## `geom_smooth()` using formula = 'y ~ x'
ggsave('pupaeinteraction2021revised.png',plot=venombyratio,dpi=900,width =30, height = 15, units = "cm",device='png')## `geom_smooth()` using formula = 'y ~ x'Modeling Worker data
by worker
Next, lets run a generalized linear mixed-effects model
workerdata$Pupae.Adult.Ratio=workerdata$Pupae.Count/workerdata$Adult.Count
workerdata$Larvae.Adult.Ratio=workerdata$Larvae.Count/workerdata$Adult.Count
workerdata$Young.Adult.Ratio=(workerdata$Larvae.Count+workerdata$Pupae.Count)/workerdata$Adult.Count
model1<-lme(Venom.Volume~Callow+Collection.Day+Behavior.Code*Adult.Count.Scaled+Behavior.Code*Pupae.Adult.Ratio+as.factor(Queen.Count>0)+Larvae.Adult.Ratio+Winged.Queen.Count+I(Webers.Length.mm^3), data=workerdata, random= ~1|Colony)
summary(model1)## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 2959.889 3038.6 -1465.945
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 0.189784 0.4673751
##
## Fixed effects: Venom.Volume ~ Callow + Collection.Day + Behavior.Code * Adult.Count.Scaled + Behavior.Code * Pupae.Adult.Ratio + as.factor(Queen.Count > 0) + Larvae.Adult.Ratio + Winged.Queen.Count + I(Webers.Length.mm^3)
## Value Std.Error DF t-value
## (Intercept) 0.8836716 0.08663683 1895 10.199723
## CallowTRUE -0.5720680 0.12467505 1895 -4.588472
## Collection.Day -0.0021061 0.00043644 148 -4.825719
## Behavior.CodeNurse -0.2667093 0.02872533 1895 -9.284812
## Adult.Count.Scaled 0.0495090 0.02514650 148 1.968821
## Pupae.Adult.Ratio -0.0235322 0.03124439 148 -0.753167
## as.factor(Queen.Count > 0)TRUE -0.0736151 0.04638146 148 -1.587167
## Larvae.Adult.Ratio -0.0065485 0.01414176 148 -0.463058
## Winged.Queen.Count 0.0129051 0.01023619 148 1.260730
## I(Webers.Length.mm^3) 1.5673781 0.14040340 1895 11.163392
## Behavior.CodeNurse:Adult.Count.Scaled -0.0451118 0.01819470 1895 -2.479390
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.1258156 0.03079669 1895 4.085362
## p-value
## (Intercept) 0.0000
## CallowTRUE 0.0000
## Collection.Day 0.0000
## Behavior.CodeNurse 0.0000
## Adult.Count.Scaled 0.0508
## Pupae.Adult.Ratio 0.4525
## as.factor(Queen.Count > 0)TRUE 0.1146
## Larvae.Adult.Ratio 0.6440
## Winged.Queen.Count 0.2094
## I(Webers.Length.mm^3) 0.0000
## Behavior.CodeNurse:Adult.Count.Scaled 0.0132
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.0000
## Correlation:
## (Intr) ClTRUE Cllc.D Bhv.CN Ad.C.S Pp.A.R
## CallowTRUE -0.026
## Collection.Day -0.353 -0.009
## Behavior.CodeNurse -0.156 -0.045 -0.036
## Adult.Count.Scaled 0.061 0.002 -0.403 0.251
## Pupae.Adult.Ratio -0.197 0.011 0.009 0.324 0.259
## as.factor(Queen.Count > 0)TRUE -0.390 0.014 -0.109 -0.024 -0.121 -0.044
## Larvae.Adult.Ratio -0.268 0.018 0.125 -0.007 0.201 0.153
## Winged.Queen.Count -0.066 -0.035 -0.050 0.019 0.074 0.005
## I(Webers.Length.mm^3) -0.735 0.029 -0.017 -0.063 -0.040 0.007
## Behavior.CodeNurse:Adult.Count.Scaled 0.093 0.006 0.013 -0.476 -0.464 -0.185
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.120 -0.015 0.025 -0.482 -0.142 -0.690
## a.(Q>0 Lr.A.R Wn.Q.C I(W.L. B.CN:A
## CallowTRUE
## Collection.Day
## Behavior.CodeNurse
## Adult.Count.Scaled
## Pupae.Adult.Ratio
## as.factor(Queen.Count > 0)TRUE
## Larvae.Adult.Ratio -0.129
## Winged.Queen.Count 0.073 0.129
## I(Webers.Length.mm^3) 0.114 -0.034 -0.033
## Behavior.CodeNurse:Adult.Count.Scaled 0.000 0.016 -0.007 0.007
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.021 0.007 0.013 -0.020 0.284
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.95840123 -0.64120107 -0.06105864 0.56990493 4.08541790
##
## Number of Observations: 2055
## Number of Groups: 155emtrends(model1,pairwise~Behavior.Code, var = "Pupae.Adult.Ratio",infer=TRUE)## $emtrends
## Behavior.Code Pupae.Adult.Ratio.trend SE df lower.CL upper.CL t.ratio
## Forager -0.0235 0.0312 148 -0.0853 0.0382 -0.753
## Nurse 0.1023 0.0244 148 0.0540 0.1505 4.190
## p.value
## 0.4525
## <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse -0.126 0.0308 1895 -0.186 -0.0654 -4.085 <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emtrends(model1,pairwise~Behavior.Code, var = "Adult.Count.Scaled",infer=TRUE)## $emtrends
## Behavior.Code Adult.Count.Scaled.trend SE df lower.CL upper.CL t.ratio
## Forager 0.0495 0.0251 148 -0.000184 0.0992 1.969
## Nurse 0.0044 0.0232 148 -0.041468 0.0503 0.189
## p.value
## 0.0508
## 0.8500
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse 0.0451 0.0182 1895 0.00943 0.0808 2.479 0.0132
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95#Remove Larvae.Adult.Ratio from model
model2<-lme(Venom.Volume~Callow+Collection.Day+Behavior.Code*Adult.Count.Scaled+Behavior.Code*Pupae.Adult.Ratio+as.factor(Queen.Count>0)+Winged.Queen.Count, data=workerdata, random= ~1|Colony)
summary(model2)## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 3066.233 3133.711 -1521.117
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 0.2199257 0.4781123
##
## Fixed effects: Venom.Volume ~ Callow + Collection.Day + Behavior.Code * Adult.Count.Scaled + Behavior.Code * Pupae.Adult.Ratio + as.factor(Queen.Count > 0) + Winged.Queen.Count
## Value Std.Error DF t-value
## (Intercept) 1.5912310 0.05814547 1896 27.366380
## CallowTRUE -0.5966223 0.12790973 1896 -4.664401
## Collection.Day -0.0020100 0.00048250 149 -4.165821
## Behavior.CodeNurse -0.2470789 0.02943246 1896 -8.394778
## Adult.Count.Scaled 0.0607093 0.02717030 149 2.234400
## Pupae.Adult.Ratio -0.0256711 0.03287725 149 -0.780817
## as.factor(Queen.Count > 0)TRUE -0.1318781 0.05075966 149 -2.598089
## Winged.Queen.Count 0.0172544 0.01138232 149 1.515894
## Behavior.CodeNurse:Adult.Count.Scaled -0.0460432 0.01865499 1896 -2.468144
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.1323307 0.03156880 1896 4.191821
## p-value
## (Intercept) 0.0000
## CallowTRUE 0.0000
## Collection.Day 0.0001
## Behavior.CodeNurse 0.0000
## Adult.Count.Scaled 0.0269
## Pupae.Adult.Ratio 0.4361
## as.factor(Queen.Count > 0)TRUE 0.0103
## Winged.Queen.Count 0.1317
## Behavior.CodeNurse:Adult.Count.Scaled 0.0137
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.0000
## Correlation:
## (Intr) ClTRUE Cllc.D Bhv.CN Ad.C.S Pp.A.R
## CallowTRUE 0.001
## Collection.Day -0.549 -0.010
## Behavior.CodeNurse -0.313 -0.043 -0.033
## Adult.Count.Scaled 0.175 -0.001 -0.447 0.240
## Pupae.Adult.Ratio -0.232 0.008 -0.012 0.318 0.235
## as.factor(Queen.Count > 0)TRUE -0.573 0.011 -0.094 -0.018 -0.099 -0.032
## Winged.Queen.Count -0.088 -0.033 -0.069 0.016 0.047 -0.015
## Behavior.CodeNurse:Adult.Count.Scaled 0.159 0.006 0.010 -0.478 -0.445 -0.183
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.163 -0.014 0.022 -0.485 -0.138 -0.673
## a.(Q>0 Wn.Q.C B.CN:A
## CallowTRUE
## Collection.Day
## Behavior.CodeNurse
## Adult.Count.Scaled
## Pupae.Adult.Ratio
## as.factor(Queen.Count > 0)TRUE
## Winged.Queen.Count 0.098
## Behavior.CodeNurse:Adult.Count.Scaled 0.003 -0.008
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.016 0.011 0.285
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.07612331 -0.65741486 -0.07120184 0.58827104 4.05847942
##
## Number of Observations: 2055
## Number of Groups: 155emtrends(model2,pairwise~Behavior.Code, var = "Pupae.Adult.Ratio",infer=TRUE)## $emtrends
## Behavior.Code Pupae.Adult.Ratio.trend SE df lower.CL upper.CL t.ratio
## Forager -0.0257 0.0329 149 -0.0906 0.0393 -0.781
## Nurse 0.1067 0.0261 149 0.0551 0.1582 4.086
## p.value
## 0.4361
## 0.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse -0.132 0.0316 1896 -0.194 -0.0704 -4.192 <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emtrends(model2,pairwise~Behavior.Code, var = "Adult.Count.Scaled",infer=TRUE)## $emtrends
## Behavior.Code Adult.Count.Scaled.trend SE df lower.CL upper.CL t.ratio
## Forager 0.0607 0.0272 149 0.00702 0.1144 2.234
## Nurse 0.0147 0.0252 149 -0.03515 0.0645 0.582
## p.value
## 0.0269
## 0.5616
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse 0.046 0.0187 1896 0.00946 0.0826 2.468 0.0137
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95#Remove Winged.Queen.Count from model
model3<-lme(Venom.Volume~Callow+Behavior.Code*Adult.Count.Scaled+Collection.Day+Behavior.Code*Pupae.Adult.Ratio+as.factor(Queen.Count>0)+I(Webers.Length.mm^3), data=workerdata, random= ~1|Colony)
summary(model3)## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 2943.851 3011.329 -1459.925
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 0.188721 0.4674991
##
## Fixed effects: Venom.Volume ~ Callow + Behavior.Code * Adult.Count.Scaled + Collection.Day + Behavior.Code * Pupae.Adult.Ratio + as.factor(Queen.Count > 0) + I(Webers.Length.mm^3)
## Value Std.Error DF t-value
## (Intercept) 0.8762991 0.08330639 1895 10.518990
## CallowTRUE -0.5654142 0.12457937 1895 -4.538586
## Behavior.CodeNurse -0.2675521 0.02872156 1895 -9.315377
## Adult.Count.Scaled 0.0502497 0.02452652 150 2.048789
## Collection.Day -0.0020420 0.00043046 150 -4.743735
## Pupae.Adult.Ratio -0.0206642 0.03082660 150 -0.670337
## as.factor(Queen.Count > 0)TRUE -0.0819851 0.04564534 150 -1.796132
## I(Webers.Length.mm^3) 1.5713074 0.14021623 1895 11.206316
## Behavior.CodeNurse:Adult.Count.Scaled -0.0447777 0.01819426 1895 -2.461087
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.1254461 0.03079856 1895 4.073115
## p-value
## (Intercept) 0.0000
## CallowTRUE 0.0000
## Behavior.CodeNurse 0.0000
## Adult.Count.Scaled 0.0422
## Collection.Day 0.0000
## Pupae.Adult.Ratio 0.5037
## as.factor(Queen.Count > 0)TRUE 0.0745
## I(Webers.Length.mm^3) 0.0000
## Behavior.CodeNurse:Adult.Count.Scaled 0.0139
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.0000
## Correlation:
## (Intr) ClTRUE Bhv.CN Ad.C.S Cllc.D Pp.A.R
## CallowTRUE -0.024
## Behavior.CodeNurse -0.163 -0.044
## Adult.Count.Scaled 0.123 0.000 0.258
## Collection.Day -0.336 -0.014 -0.034 -0.439
## Pupae.Adult.Ratio -0.164 0.008 0.329 0.237 -0.012
## as.factor(Queen.Count > 0)TRUE -0.443 0.020 -0.027 -0.103 -0.089 -0.023
## I(Webers.Length.mm^3) -0.776 0.028 -0.063 -0.033 -0.015 0.012
## Behavior.CodeNurse:Adult.Count.Scaled 0.101 0.006 -0.476 -0.479 0.010 -0.190
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.127 -0.014 -0.483 -0.148 0.025 -0.701
## a.(Q>0 I(W.L. B.CN:A
## CallowTRUE
## Behavior.CodeNurse
## Adult.Count.Scaled
## Collection.Day
## Pupae.Adult.Ratio
## as.factor(Queen.Count > 0)TRUE
## I(Webers.Length.mm^3) 0.114
## Behavior.CodeNurse:Adult.Count.Scaled 0.003 0.007
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.022 -0.019 0.284
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.96693601 -0.64399389 -0.06399286 0.57044886 4.08949506
##
## Number of Observations: 2055
## Number of Groups: 155emtrends(model3,pairwise~Behavior.Code, var = "Pupae.Adult.Ratio",infer=TRUE)## $emtrends
## Behavior.Code Pupae.Adult.Ratio.trend SE df lower.CL upper.CL t.ratio
## Forager -0.0207 0.0308 150 -0.0816 0.0402 -0.670
## Nurse 0.1048 0.0238 150 0.0577 0.1519 4.397
## p.value
## 0.5037
## <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse -0.125 0.0308 1895 -0.186 -0.065 -4.073 <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emtrends(model3,pairwise~Behavior.Code, var = "Adult.Count.Scaled",infer=TRUE)## $emtrends
## Behavior.Code Adult.Count.Scaled.trend SE df lower.CL upper.CL t.ratio
## Forager 0.05025 0.0245 150 0.00179 0.0987 2.049
## Nurse 0.00547 0.0225 150 -0.03894 0.0499 0.243
## p.value
## 0.0422
## 0.8080
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse 0.0448 0.0182 1895 0.00909 0.0805 2.461 0.0139
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95#Callow workers have significantly less venom
#Workers have less venom as season progresses from may to october
#significant interaction effect between behavior and adult count
#significant interaction effect between behavior and pupae count
plot(model3)
qqnorm(resid(model3))
hist(resid(model3))
#Overall, slopes are significantly different from each other. MEANING that the slope for nurses to predict venom based on adult count is significantly different than that of foragers, and same for pupae.
#For Pupae adult ratio, slope for nurses is the only one that is significantly different than 0. So as the number of pupae in a colony increase, the nurses in that colony have more venom.
#This is another way to do an ad hoc test below. Basically, this code says at -1,0,and 1, slopes for nurses and foragers are significantly different from each other. Here it is only written for worker count
pairs(emmeans(model3, ~ Behavior.Code * Adult.Count.Scaled,at=list(Adult.Count.Scaled=c(-1,0,1))))## contrast estimate SE
## (Forager Adult.Count.Scaled-1) - (Nurse Adult.Count.Scaled-1) 0.17511 0.0367
## (Forager Adult.Count.Scaled-1) - Forager Adult.Count.Scaled0 -0.05025 0.0245
## (Forager Adult.Count.Scaled-1) - Nurse Adult.Count.Scaled0 0.16964 0.0390
## (Forager Adult.Count.Scaled-1) - Forager Adult.Count.Scaled1 -0.10050 0.0491
## (Forager Adult.Count.Scaled-1) - Nurse Adult.Count.Scaled1 0.16417 0.0520
## (Nurse Adult.Count.Scaled-1) - Forager Adult.Count.Scaled0 -0.22536 0.0352
## (Nurse Adult.Count.Scaled-1) - Nurse Adult.Count.Scaled0 -0.00547 0.0225
## (Nurse Adult.Count.Scaled-1) - Forager Adult.Count.Scaled1 -0.27561 0.0484
## (Nurse Adult.Count.Scaled-1) - Nurse Adult.Count.Scaled1 -0.01094 0.0450
## Forager Adult.Count.Scaled0 - Nurse Adult.Count.Scaled0 0.21989 0.0252
## Forager Adult.Count.Scaled0 - Forager Adult.Count.Scaled1 -0.05025 0.0245
## Forager Adult.Count.Scaled0 - Nurse Adult.Count.Scaled1 0.21442 0.0323
## Nurse Adult.Count.Scaled0 - Forager Adult.Count.Scaled1 -0.27014 0.0310
## Nurse Adult.Count.Scaled0 - Nurse Adult.Count.Scaled1 -0.00547 0.0225
## Forager Adult.Count.Scaled1 - Nurse Adult.Count.Scaled1 0.26467 0.0243
## df t.ratio p.value
## 1895 4.776 <.0001
## 150 -2.049 0.3201
## 150 4.355 0.0003
## 150 -2.049 0.3201
## 150 3.159 0.0231
## 150 -6.399 <.0001
## 150 -0.243 0.9999
## 150 -5.698 <.0001
## 150 -0.243 0.9999
## 1895 8.710 <.0001
## 150 -2.049 0.3201
## 150 6.633 <.0001
## 150 -8.716 <.0001
## 150 -0.243 0.9999
## 1895 10.875 <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 6 estimatespairs(emmeans(model3, ~ Behavior.Code | Adult.Count.Scaled,at=list(Adult.Count.Scaled=c(-1,0,1))))## Adult.Count.Scaled = -1:
## contrast estimate SE df t.ratio p.value
## Forager - Nurse 0.175 0.0367 1895 4.776 <.0001
##
## Adult.Count.Scaled = 0:
## contrast estimate SE df t.ratio p.value
## Forager - Nurse 0.220 0.0252 1895 8.710 <.0001
##
## Adult.Count.Scaled = 1:
## contrast estimate SE df t.ratio p.value
## Forager - Nurse 0.265 0.0243 1895 10.875 <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containmentrsq.lmm(model3)## $model
## [1] 0.3104088
##
## $fixed
## [1] 0.2001303
##
## $random
## [1] 0.1102786In this chunk, I am using the performance package function check_collinearity to test if date and pupae are collinear (and a problem), since pupae are only existant at certain points during the season
check_collinearity(model3)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## Callow 1.00 [1.00, 42.71] 1.00 1.00
## Behavior.Code 1.58 [1.49, 1.68] 1.26 0.63
## Adult.Count.Scaled 1.90 [1.79, 2.03] 1.38 0.53
## Collection.Day 1.39 [1.32, 1.47] 1.18 0.72
## Pupae.Adult.Ratio 2.11 [1.98, 2.25] 1.45 0.47
## as.factor(Queen.Count > 0) 1.06 [1.02, 1.13] 1.03 0.95
## I(Webers.Length.mm^3) 1.02 [1.00, 1.17] 1.01 0.98
## Behavior.Code:Adult.Count.Scaled 1.72 [1.62, 1.83] 1.31 0.58
## Behavior.Code:Pupae.Adult.Ratio 2.38 [2.22, 2.55] 1.54 0.42
## Tolerance 95% CI
## [0.02, 1.00]
## [0.60, 0.67]
## [0.49, 0.56]
## [0.68, 0.76]
## [0.44, 0.51]
## [0.88, 0.98]
## [0.85, 1.00]
## [0.55, 0.62]
## [0.39, 0.45]#No problems here! All well under 5, and even under 3.###Check Outliers We found a few nests with very high pupae to worker ratios. If we remove these outliers, do our results stay the same?
#What the figure looks like with outliers removed (Pupae.Adult.Ratio>=5 removed)
venombyratio2=ggplot(data=subset(workerdata,Pupae.Adult.Ratio<5), aes(Pupae.Adult.Ratio,standardizedvenom, color=Behavior.Code)) +
geom_point() +
geom_smooth(method=lm, fill="gray", se=TRUE) +
xlab("Ratio of Pupae to Adults")+
ylab("Standardized Venom Volume")+
theme_bw(base_size = 22)+
scale_colour_discrete("Behavior")
venombyratio2## `geom_smooth()` using formula = 'y ~ x'
#Model results with removed
model3_OutliersRemoved<-lme(Venom.Volume~Callow+Behavior.Code*Adult.Count.Scaled+Collection.Day+Behavior.Code*Pupae.Adult.Ratio+as.factor(Queen.Count>0)+I(Webers.Length.mm^3), data=subset(workerdata,Pupae.Adult.Ratio<5), random= ~1|Colony)
summary(model3_OutliersRemoved)## Linear mixed-effects model fit by REML
## Data: subset(workerdata, Pupae.Adult.Ratio < 5)
## AIC BIC logLik
## 2934.333 3001.782 -1455.167
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 0.1894164 0.4671922
##
## Fixed effects: Venom.Volume ~ Callow + Behavior.Code * Adult.Count.Scaled + Collection.Day + Behavior.Code * Pupae.Adult.Ratio + as.factor(Queen.Count > 0) + I(Webers.Length.mm^3)
## Value Std.Error DF t-value
## (Intercept) 0.8836272 0.08388724 1891 10.533512
## CallowTRUE -0.5673943 0.12452186 1891 -4.556584
## Behavior.CodeNurse -0.2760435 0.02906635 1891 -9.497011
## Adult.Count.Scaled 0.0484096 0.02459451 149 1.968308
## Collection.Day -0.0020381 0.00043208 149 -4.716903
## Pupae.Adult.Ratio -0.0417924 0.03514675 149 -1.189084
## as.factor(Queen.Count > 0)TRUE -0.0831602 0.04576463 149 -1.817129
## I(Webers.Length.mm^3) 1.5705031 0.14026564 1891 11.196635
## Behavior.CodeNurse:Adult.Count.Scaled -0.0423318 0.01823115 1891 -2.321951
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.1527847 0.03417150 1891 4.471116
## p-value
## (Intercept) 0.0000
## CallowTRUE 0.0000
## Behavior.CodeNurse 0.0000
## Adult.Count.Scaled 0.0509
## Collection.Day 0.0000
## Pupae.Adult.Ratio 0.2363
## as.factor(Queen.Count > 0)TRUE 0.0712
## I(Webers.Length.mm^3) 0.0000
## Behavior.CodeNurse:Adult.Count.Scaled 0.0203
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.0000
## Correlation:
## (Intr) ClTRUE Bhv.CN Ad.C.S Cllc.D Pp.A.R
## CallowTRUE -0.024
## Behavior.CodeNurse -0.170 -0.042
## Adult.Count.Scaled 0.119 0.001 0.260
## Collection.Day -0.339 -0.014 -0.034 -0.437
## Pupae.Adult.Ratio -0.195 0.009 0.339 0.226 0.008
## as.factor(Queen.Count > 0)TRUE -0.444 0.020 -0.024 -0.102 -0.087 -0.005
## I(Webers.Length.mm^3) -0.773 0.028 -0.060 -0.032 -0.014 0.023
## Behavior.CodeNurse:Adult.Count.Scaled 0.105 0.005 -0.480 -0.479 0.010 -0.191
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.142 -0.018 -0.497 -0.151 0.022 -0.710
## a.(Q>0 I(W.L. B.CN:A
## CallowTRUE
## Behavior.CodeNurse
## Adult.Count.Scaled
## Collection.Day
## Pupae.Adult.Ratio
## as.factor(Queen.Count > 0)TRUE
## I(Webers.Length.mm^3) 0.115
## Behavior.CodeNurse:Adult.Count.Scaled 0.002 0.006
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.027 -0.024 0.286
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.96902240 -0.64210328 -0.05755629 0.56627334 4.08565791
##
## Number of Observations: 2050
## Number of Groups: 154emtrends(model3_OutliersRemoved,pairwise~Behavior.Code, var = "Pupae.Adult.Ratio",infer=TRUE)## $emtrends
## Behavior.Code Pupae.Adult.Ratio.trend SE df lower.CL upper.CL t.ratio
## Forager -0.0418 0.0351 149 -0.1112 0.0277 -1.189
## Nurse 0.1110 0.0264 149 0.0588 0.1632 4.202
## p.value
## 0.2363
## <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse -0.153 0.0342 1891 -0.22 -0.0858 -4.471 <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emtrends(model3_OutliersRemoved,pairwise~Behavior.Code, var = "Adult.Count.Scaled",infer=TRUE)## $emtrends
## Behavior.Code Adult.Count.Scaled.trend SE df lower.CL upper.CL t.ratio
## Forager 0.04841 0.0246 149 -0.00019 0.0970 1.968
## Nurse 0.00608 0.0225 149 -0.03845 0.0506 0.270
## p.value
## 0.0509
## 0.7877
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse 0.0423 0.0182 1891 0.00658 0.0781 2.322 0.0203
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95#Nurses from colonies with higher pupae to adult ratios still have more venom even when we remove these outliers###Controlling for Size of ants Which is the best way to control for Weber’s length? Do we standardize response variable or include it in the model? Do we cube it or leave it linear? Some of this code from Lynn Johnson at cornell statistical consulting unit
plot(workerdata$standardizedvenom, workerdata$Venom.Volume/workerdata$Webers.Length.mm)
workerdata$standardizedvenom2 <- workerdata$Venom.Volume/(workerdata$Webers.Length.mm^3)
workerdata$Pupae.Adult.Ratio <- workerdata$Pupae.Count/workerdata$Adult.Count
# venom/weber
model3 <- lme(standardizedvenom ~ Callow+Behavior.Code*Adult.Count.Scaled+Collection.Day+
Behavior.Code*Pupae.Adult.Ratio+as.factor(Queen.Count>0), random= ~1|Colony,
data=workerdata)
summary(model3)## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 4073.704 4135.564 -2025.852
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 0.2624909 0.6150299
##
## Fixed effects: standardizedvenom ~ Callow + Behavior.Code * Adult.Count.Scaled + Collection.Day + Behavior.Code * Pupae.Adult.Ratio + as.factor(Queen.Count > 0)
## Value Std.Error DF t-value
## (Intercept) 2.0838899 0.07132631 1896 29.216285
## CallowTRUE -0.7781286 0.16410284 1896 -4.741713
## Behavior.CodeNurse -0.3417640 0.03777374 1896 -9.047663
## Adult.Count.Scaled 0.0708162 0.03333070 150 2.124655
## Collection.Day -0.0026611 0.00058777 150 -4.527375
## Pupae.Adult.Ratio -0.0350806 0.04124894 150 -0.850460
## as.factor(Queen.Count > 0)TRUE -0.1411956 0.06182260 150 -2.283883
## Behavior.CodeNurse:Adult.Count.Scaled -0.0591867 0.02396180 1896 -2.470043
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.1727077 0.04055188 1896 4.258932
## p-value
## (Intercept) 0.0000
## CallowTRUE 0.0000
## Behavior.CodeNurse 0.0000
## Adult.Count.Scaled 0.0353
## Collection.Day 0.0000
## Pupae.Adult.Ratio 0.3964
## as.factor(Queen.Count > 0)TRUE 0.0238
## Behavior.CodeNurse:Adult.Count.Scaled 0.0136
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.0000
## Correlation:
## (Intr) ClTRUE Bhv.CN Ad.C.S Cllc.D Pp.A.R
## CallowTRUE -0.002
## Behavior.CodeNurse -0.327 -0.043
## Adult.Count.Scaled 0.166 0.001 0.249
## Collection.Day -0.555 -0.013 -0.034 -0.442
## Pupae.Adult.Ratio -0.241 0.008 0.326 0.237 -0.012
## as.factor(Queen.Count > 0)TRUE -0.567 0.016 -0.020 -0.102 -0.088 -0.027
## Behavior.CodeNurse:Adult.Count.Scaled 0.165 0.006 -0.477 -0.464 0.010 -0.187
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.172 -0.014 -0.485 -0.145 0.024 -0.689
## a.(Q>0 B.CN:A
## CallowTRUE
## Behavior.CodeNurse
## Adult.Count.Scaled
## Collection.Day
## Pupae.Adult.Ratio
## as.factor(Queen.Count > 0)TRUE
## Behavior.CodeNurse:Adult.Count.Scaled 0.003
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.019 0.285
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.99671084 -0.65095358 -0.07172188 0.57447689 4.02938992
##
## Number of Observations: 2055
## Number of Groups: 155emtrends(model3, pairwise~Behavior.Code, var = "Pupae.Adult.Ratio",infer=TRUE)## $emtrends
## Behavior.Code Pupae.Adult.Ratio.trend SE df lower.CL upper.CL t.ratio
## Forager -0.0351 0.0412 150 -0.1166 0.0464 -0.850
## Nurse 0.1376 0.0323 150 0.0739 0.2014 4.267
## p.value
## 0.3964
## <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse -0.173 0.0406 1896 -0.252 -0.0932 -4.259 <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emtrends(model3, pairwise~Behavior.Code, var = "Adult.Count.Scaled",infer=TRUE)## $emtrends
## Behavior.Code Adult.Count.Scaled.trend SE df lower.CL upper.CL t.ratio
## Forager 0.0708 0.0333 150 0.00496 0.1367 2.125
## Nurse 0.0116 0.0307 150 -0.04906 0.0723 0.379
## p.value
## 0.0353
## 0.7055
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse 0.0592 0.024 1896 0.0122 0.106 2.470 0.0136
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95# venom ~ weber
model3_fe <- lme(Venom.Volume ~ Callow+Behavior.Code*Adult.Count.Scaled+Collection.Day+
Behavior.Code*Pupae.Adult.Ratio+as.factor(Queen.Count>0) + Webers.Length.mm,
random= ~1|Colony, data=workerdata)
summary(model3_fe)## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 2941.951 3009.428 -1458.975
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 0.1890786 0.4673611
##
## Fixed effects: Venom.Volume ~ Callow + Behavior.Code * Adult.Count.Scaled + Collection.Day + Behavior.Code * Pupae.Adult.Ratio + as.factor(Queen.Count > 0) + Webers.Length.mm
## Value Std.Error DF t-value
## (Intercept) -0.5038156 0.19445914 1895 -2.590856
## CallowTRUE -0.5652206 0.12455057 1895 -4.538081
## Behavior.CodeNurse -0.2666903 0.02870993 1895 -9.289133
## Adult.Count.Scaled 0.0492019 0.02455356 150 2.003861
## Collection.Day -0.0020307 0.00043094 150 -4.712377
## Pupae.Adult.Ratio -0.0209936 0.03083711 150 -0.680789
## as.factor(Queen.Count > 0)TRUE -0.0824510 0.04568753 150 -1.804672
## Webers.Length.mm 2.7351602 0.24332905 1895 11.240582
## Behavior.CodeNurse:Adult.Count.Scaled -0.0454242 0.01818935 1895 -2.497298
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.1246355 0.03079217 1895 4.047636
## p-value
## (Intercept) 0.0096
## CallowTRUE 0.0000
## Behavior.CodeNurse 0.0000
## Adult.Count.Scaled 0.0469
## Collection.Day 0.0000
## Pupae.Adult.Ratio 0.4971
## as.factor(Queen.Count > 0)TRUE 0.0731
## Webers.Length.mm 0.0000
## Behavior.CodeNurse:Adult.Count.Scaled 0.0126
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.0001
## Correlation:
## (Intr) ClTRUE Bhv.CN Ad.C.S Cllc.D Pp.A.R
## CallowTRUE -0.028
## Behavior.CodeNurse -0.033 -0.044
## Adult.Count.Scaled 0.077 0.000 0.258
## Collection.Day -0.137 -0.014 -0.034 -0.439
## Pupae.Adult.Ratio -0.077 0.008 0.329 0.237 -0.012
## as.factor(Queen.Count > 0)TRUE -0.261 0.020 -0.027 -0.103 -0.089 -0.023
## Webers.Length.mm -0.963 0.028 -0.060 -0.036 -0.013 0.011
## Behavior.CodeNurse:Adult.Count.Scaled 0.042 0.006 -0.476 -0.478 0.010 -0.190
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.069 -0.015 -0.483 -0.148 0.025 -0.701
## a.(Q>0 Wbr.L. B.CN:A
## CallowTRUE
## Behavior.CodeNurse
## Adult.Count.Scaled
## Collection.Day
## Pupae.Adult.Ratio
## as.factor(Queen.Count > 0)TRUE
## Webers.Length.mm 0.113
## Behavior.CodeNurse:Adult.Count.Scaled 0.003 0.004
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.022 -0.022 0.284
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.88441720 -0.64464959 -0.06859058 0.57603452 4.06452206
##
## Number of Observations: 2055
## Number of Groups: 155emtrends(model3_fe, pairwise~Behavior.Code, var = "Pupae.Adult.Ratio",infer=TRUE)## $emtrends
## Behavior.Code Pupae.Adult.Ratio.trend SE df lower.CL upper.CL t.ratio
## Forager -0.021 0.0308 150 -0.0819 0.0399 -0.681
## Nurse 0.104 0.0238 150 0.0565 0.1508 4.346
## p.value
## 0.4971
## <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse -0.125 0.0308 1895 -0.185 -0.0642 -4.048 0.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emtrends(model3_fe, pairwise~Behavior.Code, var = "Adult.Count.Scaled",infer=TRUE)## $emtrends
## Behavior.Code Adult.Count.Scaled.trend SE df lower.CL upper.CL t.ratio
## Forager 0.04920 0.0246 150 0.000686 0.0977 2.004
## Nurse 0.00378 0.0225 150 -0.040697 0.0483 0.168
## p.value
## 0.0469
## 0.8669
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse 0.0454 0.0182 1895 0.00975 0.0811 2.497 0.0126
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95# venom/weber^3
model3_cubed <- lme(standardizedvenom2 ~ Callow+Behavior.Code*Adult.Count.Scaled+Collection.Day+
Behavior.Code*Pupae.Adult.Ratio+as.factor(Queen.Count>0), random= ~1|Colony, data=workerdata)
summary(model3_cubed)## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 6452.593 6514.453 -3215.296
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 0.4370534 1.103779
##
## Fixed effects: standardizedvenom2 ~ Callow + Behavior.Code * Adult.Count.Scaled + Collection.Day + Behavior.Code * Pupae.Adult.Ratio + as.factor(Queen.Count > 0)
## Value Std.Error DF t-value
## (Intercept) 3.583215 0.12283591 1896 29.170744
## CallowTRUE -1.341644 0.29384015 1896 -4.565899
## Behavior.CodeNurse -0.654052 0.06764028 1896 -9.669558
## Adult.Count.Scaled 0.095820 0.05723594 150 1.674125
## Collection.Day -0.004898 0.00100341 150 -4.881117
## Pupae.Adult.Ratio -0.071500 0.07236804 150 -0.988000
## as.factor(Queen.Count > 0)TRUE -0.099888 0.10577446 150 -0.944352
## Behavior.CodeNurse:Adult.Count.Scaled -0.099777 0.04293916 1896 -2.323691
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.297955 0.07267617 1896 4.099759
## p-value
## (Intercept) 0.0000
## CallowTRUE 0.0000
## Behavior.CodeNurse 0.0000
## Adult.Count.Scaled 0.0962
## Collection.Day 0.0000
## Pupae.Adult.Ratio 0.3247
## as.factor(Queen.Count > 0)TRUE 0.3465
## Behavior.CodeNurse:Adult.Count.Scaled 0.0202
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.0000
## Correlation:
## (Intr) ClTRUE Bhv.CN Ad.C.S Cllc.D Pp.A.R
## CallowTRUE -0.003
## Behavior.CodeNurse -0.340 -0.042
## Adult.Count.Scaled 0.151 0.001 0.259
## Collection.Day -0.551 -0.013 -0.035 -0.439
## Pupae.Adult.Ratio -0.247 0.008 0.333 0.237 -0.011
## as.factor(Queen.Count > 0)TRUE -0.565 0.017 -0.020 -0.099 -0.088 -0.023
## Behavior.CodeNurse:Adult.Count.Scaled 0.172 0.006 -0.476 -0.484 0.011 -0.191
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.179 -0.014 -0.485 -0.150 0.025 -0.705
## a.(Q>0 B.CN:A
## CallowTRUE
## Behavior.CodeNurse
## Adult.Count.Scaled
## Collection.Day
## Pupae.Adult.Ratio
## as.factor(Queen.Count > 0)TRUE
## Behavior.CodeNurse:Adult.Count.Scaled 0.002
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.020 0.284
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.64593742 -0.63690526 -0.05356094 0.56198871 6.07345845
##
## Number of Observations: 2055
## Number of Groups: 155emtrends(model3_cubed, pairwise~Behavior.Code, var = "Pupae.Adult.Ratio",infer=TRUE)## $emtrends
## Behavior.Code Pupae.Adult.Ratio.trend SE df lower.CL upper.CL t.ratio
## Forager -0.0715 0.0724 150 -0.214 0.0715 -0.988
## Nurse 0.2265 0.0557 150 0.116 0.3366 4.064
## p.value
## 0.3247
## 0.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse -0.298 0.0727 1896 -0.44 -0.155 -4.100 <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emtrends(model3_cubed, pairwise~Behavior.Code, var = "Adult.Count.Scaled",infer=TRUE)## $emtrends
## Behavior.Code Adult.Count.Scaled.trend SE df lower.CL upper.CL t.ratio
## Forager 0.09582 0.0572 150 -0.0173 0.2089 1.674
## Nurse -0.00396 0.0524 150 -0.1074 0.0995 -0.076
## p.value
## 0.0962
## 0.9399
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse 0.0998 0.0429 1896 0.0156 0.184 2.324 0.0202
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95# venom ~ weber^3
model3_cubed_fe <- lme(Venom.Volume ~ Callow+Behavior.Code*Adult.Count.Scaled+Collection.Day+
Behavior.Code*Pupae.Adult.Ratio+as.factor(Queen.Count>0) + I(Webers.Length.mm^3),
random= ~1|Colony, data=workerdata)
summary(model3_cubed_fe)## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 2943.851 3011.329 -1459.925
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 0.188721 0.4674991
##
## Fixed effects: Venom.Volume ~ Callow + Behavior.Code * Adult.Count.Scaled + Collection.Day + Behavior.Code * Pupae.Adult.Ratio + as.factor(Queen.Count > 0) + I(Webers.Length.mm^3)
## Value Std.Error DF t-value
## (Intercept) 0.8762991 0.08330639 1895 10.518990
## CallowTRUE -0.5654142 0.12457937 1895 -4.538586
## Behavior.CodeNurse -0.2675521 0.02872156 1895 -9.315377
## Adult.Count.Scaled 0.0502497 0.02452652 150 2.048789
## Collection.Day -0.0020420 0.00043046 150 -4.743735
## Pupae.Adult.Ratio -0.0206642 0.03082660 150 -0.670337
## as.factor(Queen.Count > 0)TRUE -0.0819851 0.04564534 150 -1.796132
## I(Webers.Length.mm^3) 1.5713074 0.14021623 1895 11.206316
## Behavior.CodeNurse:Adult.Count.Scaled -0.0447777 0.01819426 1895 -2.461087
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.1254461 0.03079856 1895 4.073115
## p-value
## (Intercept) 0.0000
## CallowTRUE 0.0000
## Behavior.CodeNurse 0.0000
## Adult.Count.Scaled 0.0422
## Collection.Day 0.0000
## Pupae.Adult.Ratio 0.5037
## as.factor(Queen.Count > 0)TRUE 0.0745
## I(Webers.Length.mm^3) 0.0000
## Behavior.CodeNurse:Adult.Count.Scaled 0.0139
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.0000
## Correlation:
## (Intr) ClTRUE Bhv.CN Ad.C.S Cllc.D Pp.A.R
## CallowTRUE -0.024
## Behavior.CodeNurse -0.163 -0.044
## Adult.Count.Scaled 0.123 0.000 0.258
## Collection.Day -0.336 -0.014 -0.034 -0.439
## Pupae.Adult.Ratio -0.164 0.008 0.329 0.237 -0.012
## as.factor(Queen.Count > 0)TRUE -0.443 0.020 -0.027 -0.103 -0.089 -0.023
## I(Webers.Length.mm^3) -0.776 0.028 -0.063 -0.033 -0.015 0.012
## Behavior.CodeNurse:Adult.Count.Scaled 0.101 0.006 -0.476 -0.479 0.010 -0.190
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.127 -0.014 -0.483 -0.148 0.025 -0.701
## a.(Q>0 I(W.L. B.CN:A
## CallowTRUE
## Behavior.CodeNurse
## Adult.Count.Scaled
## Collection.Day
## Pupae.Adult.Ratio
## as.factor(Queen.Count > 0)TRUE
## I(Webers.Length.mm^3) 0.114
## Behavior.CodeNurse:Adult.Count.Scaled 0.003 0.007
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.022 -0.019 0.284
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.96693601 -0.64399389 -0.06399286 0.57044886 4.08949506
##
## Number of Observations: 2055
## Number of Groups: 155emtrends(model3_cubed_fe, pairwise~Behavior.Code, var = "Pupae.Adult.Ratio",infer=TRUE)## $emtrends
## Behavior.Code Pupae.Adult.Ratio.trend SE df lower.CL upper.CL t.ratio
## Forager -0.0207 0.0308 150 -0.0816 0.0402 -0.670
## Nurse 0.1048 0.0238 150 0.0577 0.1519 4.397
## p.value
## 0.5037
## <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse -0.125 0.0308 1895 -0.186 -0.065 -4.073 <.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emtrends(model3_cubed_fe, pairwise~Behavior.Code, var = "Adult.Count.Scaled",infer=TRUE)## $emtrends
## Behavior.Code Adult.Count.Scaled.trend SE df lower.CL upper.CL t.ratio
## Forager 0.05025 0.0245 150 0.00179 0.0987 2.049
## Nurse 0.00547 0.0225 150 -0.03894 0.0499 0.243
## p.value
## 0.0422
## 0.8080
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse 0.0448 0.0182 1895 0.00909 0.0805 2.461 0.0139
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95plot(workerdata$standardizedvenom, fitted(model3), xlim = c(0,0.0035), ylim = c(0,0.0035))
cor(workerdata$standardizedvenom, fitted(model3))## [1] 0.5463299plot(workerdata$Venom.Volume, fitted(model3_fe), xlim = c(0,0.0035), ylim = c(0,0.0035))
cor(workerdata$Venom.Volume, fitted(model3_fe))## [1] 0.5820932plot(workerdata$standardizedvenom2, fitted(model3_cubed), xlim = c(0,0.01), ylim = c(0,0.01))
cor(workerdata$standardizedvenom2, fitted(model3_cubed))## [1] 0.5272472plot(workerdata$Venom.Volume, fitted(model3_cubed_fe), xlim = c(0,0.004), ylim = c(0,0.004))
cor(workerdata$Venom.Volume, fitted(model3_cubed_fe))## [1] 0.5817143Breakpoint
Investigating the idea of a breakpoint analysis. Ultimately decided to use date as a continuous variable because it’s hard to determine how to bin date.
#Pupae first in colonies on june 14. Almost none after august 20
for (row in 1:nrow(workerdata)) {
if(workerdata$Date.Collected.D.M[row]<"2021-06-14") {
workerdata$Raid.Potential[row]="too early"
}else if(workerdata$Date.Collected.D.M[row]>"2021-08-20"){workerdata$Raid.Potential[row]="too late"
}else{workerdata$Raid.Potential[row]="raids possible"}
}
workerdata$Raid.Potential <- factor(workerdata$Raid.Potential, levels = c("too early", "raids possible", "too late"))
breakpointmodel<-lme(standardizedvenom~Callow+Raid.Potential+Behavior.Code*Adult.Count.Scaled+Behavior.Code*Pupae.Adult.Ratio+as.factor(Queen.Count>0), data=workerdata, random= ~1|Colony)
summary(breakpointmodel)## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 4072.962 4140.44 -2024.481
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 0.2683308 0.6150688
##
## Fixed effects: standardizedvenom ~ Callow + Raid.Potential + Behavior.Code * Adult.Count.Scaled + Behavior.Code * Pupae.Adult.Ratio + as.factor(Queen.Count > 0)
## Value Std.Error DF t-value
## (Intercept) 2.0240813 0.07303700 1896 27.713093
## CallowTRUE -0.7923336 0.16440064 1896 -4.819529
## Raid.Potentialraids possible -0.0867576 0.07881765 149 -1.100738
## Raid.Potentialtoo late -0.3184087 0.08121220 149 -3.920700
## Behavior.CodeNurse -0.3404406 0.03783028 1896 -8.999155
## Adult.Count.Scaled 0.0670000 0.03405498 149 1.967405
## Pupae.Adult.Ratio -0.0474365 0.04639849 149 -1.022371
## as.factor(Queen.Count > 0)TRUE -0.1456756 0.06267870 149 -2.324164
## Behavior.CodeNurse:Adult.Count.Scaled -0.0592683 0.02397833 1896 -2.471745
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.1730084 0.04057373 1896 4.264050
## p-value
## (Intercept) 0.0000
## CallowTRUE 0.0000
## Raid.Potentialraids possible 0.2728
## Raid.Potentialtoo late 0.0001
## Behavior.CodeNurse 0.0000
## Adult.Count.Scaled 0.0510
## Pupae.Adult.Ratio 0.3083
## as.factor(Queen.Count > 0)TRUE 0.0215
## Behavior.CodeNurse:Adult.Count.Scaled 0.0135
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.0000
## Correlation:
## (Intr) ClTRUE Rd.Ptp Rd.Ptl Bhv.CN Ad.C.S
## CallowTRUE 0.005
## Raid.Potentialraids possible -0.496 -0.042
## Raid.Potentialtoo late -0.491 0.000 0.503
## Behavior.CodeNurse -0.309 -0.042 -0.036 -0.050
## Adult.Count.Scaled 0.101 -0.011 -0.104 -0.437 0.250
## Pupae.Adult.Ratio -0.061 0.028 -0.411 -0.059 0.298 0.171
## as.factor(Queen.Count > 0)TRUE -0.585 0.013 -0.007 -0.072 -0.020 -0.104
## Behavior.CodeNurse:Adult.Count.Scaled 0.154 0.005 0.021 0.017 -0.477 -0.456
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.166 -0.014 0.018 0.026 -0.485 -0.142
## Pp.A.R a.(Q>0 B.CN:A
## CallowTRUE
## Raid.Potentialraids possible
## Raid.Potentialtoo late
## Behavior.CodeNurse
## Adult.Count.Scaled
## Pupae.Adult.Ratio
## as.factor(Queen.Count > 0)TRUE -0.037
## Behavior.CodeNurse:Adult.Count.Scaled -0.174 0.003
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.617 -0.018 0.285
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.00245342 -0.64496064 -0.06928748 0.57093296 4.01715710
##
## Number of Observations: 2055
## Number of Groups: 155emtrends(breakpointmodel,pairwise~Behavior.Code, var = "Adult.Count.Scaled",infer=TRUE)## $emtrends
## Behavior.Code Adult.Count.Scaled.trend SE df lower.CL upper.CL t.ratio
## Forager 0.06700 0.0341 149 -0.000293 0.1343 1.967
## Nurse 0.00773 0.0315 149 -0.054445 0.0699 0.246
## p.value
## 0.0510
## 0.8062
##
## Results are averaged over the levels of: Callow, Raid.Potential, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse 0.0593 0.024 1896 0.0122 0.106 2.472 0.0135
##
## Results are averaged over the levels of: Callow, Raid.Potential, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emtrends(breakpointmodel,pairwise~Behavior.Code, var = "Pupae.Adult.Ratio",infer=TRUE)## $emtrends
## Behavior.Code Pupae.Adult.Ratio.trend SE df lower.CL upper.CL t.ratio
## Forager -0.0474 0.0464 149 -0.1391 0.0442 -1.022
## Nurse 0.1256 0.0384 149 0.0496 0.2015 3.267
## p.value
## 0.3083
## 0.0014
##
## Results are averaged over the levels of: Callow, Raid.Potential, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse -0.173 0.0406 1896 -0.253 -0.0934 -4.264 <.0001
##
## Results are averaged over the levels of: Callow, Raid.Potential, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emmeans(breakpointmodel, specs = pairwise ~ Raid.Potential)## $emmeans
## Raid.Potential emmean SE df lower.CL upper.CL
## too early 1.36 0.1038 149 1.15 1.56
## raids possible 1.27 0.0957 149 1.08 1.46
## too late 1.04 0.0967 149 0.85 1.23
##
## Results are averaged over the levels of: Callow, Behavior.Code, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## too early - raids possible 0.0868 0.0788 149 1.101 0.5151
## too early - too late 0.3184 0.0812 149 3.921 0.0004
## raids possible - too late 0.2317 0.0798 149 2.904 0.0118
##
## Results are averaged over the levels of: Callow, Behavior.Code, Queen.Count
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 3 estimatesvenombyraidpotential<-ggplot(workerdata,aes(Raid.Potential,standardizedvenom),color=Behavior.Code)+
geom_boxplot()+
xlab("Raid Potential")+ylab("Standardized Venom Volume")+theme_bw(base_size = 14)+
theme(legend.position="none")+
scale_colour_discrete("")
venombyraidpotential
workerdata$Raid.Potential <- factor(workerdata$Raid.Potential, levels = c("early", "maximum potential", "too late"))
for (row in 1:nrow(workerdata)) {
if(workerdata$Date.Collected.D.M[row] <"2021-06-15") {
workerdata$Raid.Potential[row]="early"
}
else if (workerdata$Date.Collected.D.M[row]>"2021-07-20"){workerdata$Raid.Potential[row]="too late"}
else {workerdata$Raid.Potential[row]="maximum potential"}
}
maxraidpotentialmodel<-lme(standardizedvenom~Callow+Raid.Potential*Behavior.Code+Behavior.Code*Adult.Count.Scaled+Behavior.Code*Pupae.Adult.Ratio+as.factor(Queen.Count>0), data=workerdata, random= ~1|Colony)
summary(maxraidpotentialmodel)## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 4071.646 4150.356 -2021.823
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 0.2550236 0.6151882
##
## Fixed effects: standardizedvenom ~ Callow + Raid.Potential * Behavior.Code + Behavior.Code * Adult.Count.Scaled + Behavior.Code * Pupae.Adult.Ratio + as.factor(Queen.Count > 0)
## Value Std.Error DF
## (Intercept) 1.9944916 0.07690699 1894
## CallowTRUE -0.7500305 0.16426972 1894
## Raid.Potentialmaximum potential -0.0529405 0.12027876 149
## Raid.Potentialtoo late -0.2419069 0.08987558 149
## Behavior.CodeNurse -0.3131962 0.06511416 1894
## Adult.Count.Scaled 0.0511010 0.03296071 149
## Pupae.Adult.Ratio -0.0424914 0.05384995 149
## as.factor(Queen.Count > 0)TRUE -0.1271376 0.06095606 149
## Raid.Potentialmaximum potential:Behavior.CodeNurse 0.1172623 0.11658552 1894
## Raid.Potentialtoo late:Behavior.CodeNurse -0.0841981 0.08587915 1894
## Behavior.CodeNurse:Adult.Count.Scaled -0.0407192 0.02714028 1894
## Behavior.CodeNurse:Pupae.Adult.Ratio 0.1171759 0.05396654 1894
## t-value p-value
## (Intercept) 25.933813 0.0000
## CallowTRUE -4.565847 0.0000
## Raid.Potentialmaximum potential -0.440149 0.6605
## Raid.Potentialtoo late -2.691576 0.0079
## Behavior.CodeNurse -4.809955 0.0000
## Adult.Count.Scaled 1.550361 0.1232
## Pupae.Adult.Ratio -0.789070 0.4313
## as.factor(Queen.Count > 0)TRUE -2.085725 0.0387
## Raid.Potentialmaximum potential:Behavior.CodeNurse 1.005805 0.3146
## Raid.Potentialtoo late:Behavior.CodeNurse -0.980425 0.3270
## Behavior.CodeNurse:Adult.Count.Scaled -1.500323 0.1337
## Behavior.CodeNurse:Pupae.Adult.Ratio 2.171269 0.0300
## Correlation:
## (Intr) ClTRUE Rd.Ptp Rd.Ptl
## CallowTRUE -0.010
## Raid.Potentialmaximum potential -0.460 -0.005
## Raid.Potentialtoo late -0.619 -0.004 0.486
## Behavior.CodeNurse -0.565 0.004 0.364 0.506
## Adult.Count.Scaled 0.183 0.000 -0.155 -0.444
## Pupae.Adult.Ratio 0.024 0.010 -0.621 -0.134
## as.factor(Queen.Count > 0)TRUE -0.529 0.020 0.009 -0.075
## Raid.Potentialmaximum potential:Behavior.CodeNurse 0.326 -0.004 -0.669 -0.338
## Raid.Potentialtoo late:Behavior.CodeNurse 0.459 -0.039 -0.342 -0.640
## Behavior.CodeNurse:Adult.Count.Scaled -0.061 0.025 0.108 0.290
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.002 -0.018 0.418 0.098
## Bhv.CN Ad.C.S Pp.A.R a.(Q>0
## CallowTRUE
## Raid.Potentialmaximum potential
## Raid.Potentialtoo late
## Behavior.CodeNurse
## Adult.Count.Scaled -0.049
## Pupae.Adult.Ratio -0.006 0.191
## as.factor(Queen.Count > 0)TRUE 0.004 -0.097 -0.040
## Raid.Potentialmaximum potential:Behavior.CodeNurse -0.550 0.093 0.434 -0.015
## Raid.Potentialtoo late:Behavior.CodeNurse -0.798 0.246 0.098 -0.020
## Behavior.CodeNurse:Adult.Count.Scaled 0.110 -0.525 -0.133 0.012
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.010 -0.113 -0.699 -0.009
## R.Pp:B R.Pl:B B.CN:A
## CallowTRUE
## Raid.Potentialmaximum potential
## Raid.Potentialtoo late
## Behavior.CodeNurse
## Adult.Count.Scaled
## Pupae.Adult.Ratio
## as.factor(Queen.Count > 0)TRUE
## Raid.Potentialmaximum potential:Behavior.CodeNurse
## Raid.Potentialtoo late:Behavior.CodeNurse 0.512
## Behavior.CodeNurse:Adult.Count.Scaled -0.161 -0.462
## Behavior.CodeNurse:Pupae.Adult.Ratio -0.618 -0.118 0.185
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.99260496 -0.64299782 -0.07551255 0.58436320 3.97873499
##
## Number of Observations: 2055
## Number of Groups: 155emtrends(maxraidpotentialmodel,pairwise~Behavior.Code, var = "Adult.Count.Scaled",infer=TRUE)## $emtrends
## Behavior.Code Adult.Count.Scaled.trend SE df lower.CL upper.CL t.ratio
## Forager 0.0511 0.0330 149 -0.0140 0.1162 1.550
## Nurse 0.0104 0.0297 149 -0.0484 0.0691 0.349
## p.value
## 0.1232
## 0.7275
##
## Results are averaged over the levels of: Callow, Raid.Potential, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse 0.0407 0.0271 1894 -0.0125 0.0939 1.500 0.1337
##
## Results are averaged over the levels of: Callow, Raid.Potential, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emtrends(maxraidpotentialmodel,pairwise~Behavior.Code, var = "Pupae.Adult.Ratio",infer=TRUE)## $emtrends
## Behavior.Code Pupae.Adult.Ratio.trend SE df lower.CL upper.CL t.ratio
## Forager -0.0425 0.0538 149 -0.14890 0.0639 -0.789
## Nurse 0.0747 0.0418 149 -0.00799 0.1574 1.785
## p.value
## 0.4313
## 0.0763
##
## Results are averaged over the levels of: Callow, Raid.Potential, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Forager - Nurse -0.117 0.054 1894 -0.223 -0.0113 -2.171 0.0300
##
## Results are averaged over the levels of: Callow, Raid.Potential, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95emmeans(maxraidpotentialmodel, specs = pairwise ~ Raid.Potential|Behavior.Code)## $emmeans
## Behavior.Code = Forager:
## Raid.Potential emmean SE df lower.CL upper.CL
## early 1.514 0.1115 149 1.294 1.73
## maximum potential 1.461 0.1214 149 1.221 1.70
## too late 1.272 0.0964 149 1.082 1.46
##
## Behavior.Code = Nurse:
## Raid.Potential emmean SE df lower.CL upper.CL
## early 1.224 0.1035 149 1.020 1.43
## maximum potential 1.289 0.1102 149 1.071 1.51
## too late 0.898 0.0896 149 0.721 1.08
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## Behavior.Code = Forager:
## contrast estimate SE df t.ratio p.value
## early - maximum potential 0.0529 0.1203 149 0.440 0.8988
## early - too late 0.2419 0.0899 149 2.692 0.0215
## maximum potential - too late 0.1890 0.1097 149 1.722 0.2003
##
## Behavior.Code = Nurse:
## contrast estimate SE df t.ratio p.value
## early - maximum potential -0.0643 0.0964 149 -0.667 0.7827
## early - too late 0.3261 0.0747 149 4.366 0.0001
## maximum potential - too late 0.3904 0.0908 149 4.299 0.0001
##
## Results are averaged over the levels of: Callow, Queen.Count
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 3 estimatesNest data
#Since dataset doesn't have year in same column, we have to jump through some hoops to get this all set up
par(mfrow=c(2,1))
plot(summarydata$Date.Collected,summarydata$Pupae.Count, xlab="Date of Colony Collection", ylab="Pupae Count")
g = ggplot(data=summarydata, aes(x=Date.Collected)) + geom_point(aes(y=Adult.Count),color="red") + geom_point(aes(y=Pupae.Count),color="blue")+ geom_point(aes(y=Larvae.Count),color="orange")+
xlab("Date Collected")+ylab("Nest Census Count")+
theme_bw(base_size = 22)
summarydata$bin <- cut(summarydata$Date.Collected, breaks=c("months"))
ggplot(data=summarydata) + geom_boxplot(aes(bin, Adult.Count),color="red")+ geom_boxplot(aes(bin, Pupae.Count),color="blue")
dat <- reshape2::melt(summarydata, id.vars="Date.Collected",
measure.vars=c("Adult.Count", "Larvae.Count","Pupae.Count"))
dat$bin <- cut(dat$Date.Collected, breaks=c("months"))
g<-ggplot(data=dat) + geom_boxplot(aes(bin, value, fill=variable))+
xlab("Date Collected")+
ylab("Nest Census")+
theme_bw(base_size = 20)+
scale_x_discrete(labels=c("May","June","July","Aug","Sep","Oct"))+
theme(legend.position = c(0.095, 0.85),legend.background = element_rect(fill = "white"))+
scale_fill_brewer(palette="BuPu",name = "Ant Type", labels = c("Workers", "Larva", "Pupa"))+
theme(legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'))
queen=ggplot(data=summarydata, aes(x=Date.Collected, y=Queen.Count)) + geom_point() +
xlab("Date Collected")+ylab("Number of Dealate Queens")+
theme_bw(base_size = 22)+
scale_x_date(date_breaks = "months" , date_labels = "%b")+
geom_smooth(method='lm', formula= y~x)
summarydata$Collection.Day<-as.numeric(summarydata$Date.Collected)-18750
censusadult<-lme(Adult.Count~bin, data=summarydata, random= ~1|Colony)
summary(censusadult)## Linear mixed-effects model fit by REML
## Data: summarydata
## AIC BIC logLik
## 1615.255 1640.151 -799.6277
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 26.38513 9.894424
##
## Fixed effects: Adult.Count ~ bin
## Value Std.Error DF t-value p-value
## (Intercept) 21.51111 4.200727 166 5.120807 0.0000
## bin2021-06-01 -3.72733 6.253606 166 -0.596029 0.5520
## bin2021-07-01 -1.02835 6.710290 166 -0.153250 0.8784
## bin2021-08-01 25.10794 7.447092 166 3.371509 0.0009
## bin2021-09-01 29.64889 7.029160 166 4.217985 0.0000
## bin2021-10-01 53.15556 8.401453 166 6.326948 0.0000
## Correlation:
## (Intr) b2021-06 b2021-07 b2021-08 b2021-09
## bin2021-06-01 -0.672
## bin2021-07-01 -0.626 0.421
## bin2021-08-01 -0.564 0.379 0.353
## bin2021-09-01 -0.598 0.401 0.374 0.337
## bin2021-10-01 -0.500 0.336 0.313 0.282 0.299
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.71854503 -0.17392770 -0.05620988 0.10237667 1.86705412
##
## Number of Observations: 172
## Number of Groups: 172#Adult count significantly increases with date, aug/sep/oct significantly higher than may/june/july
censusadult<-lme(Adult.Count~Collection.Day, data=summarydata, random= ~1|Colony)
summary(censusadult)## Linear mixed-effects model fit by REML
## Data: summarydata
## AIC BIC logLik
## 1653.608 1666.152 -822.8042
##
## Random effects:
## Formula: ~1 | Colony
## (Intercept) Residual
## StdDev: 27.16369 10.18639
##
## Fixed effects: Adult.Count ~ Collection.Day
## Value Std.Error DF t-value p-value
## (Intercept) 9.747790 3.796915 170 2.567292 0.0111
## Collection.Day 0.319956 0.043309 170 7.387799 0.0000
## Correlation:
## (Intr)
## Collection.Day -0.813
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.61558692 -0.21182705 -0.07487622 0.11299928 1.73775995
##
## Number of Observations: 172
## Number of Groups: 172